#create a vector of length 3
n0<-c(200,20,16)
n0

#create a matrix and fill it with values from 
#the vector
#c(0,.2,0,4,0,.5,2.5,0,0)
#such that the matrix should have 3 rows
M<-matrix(c(0,.2,0,4,0,.5,2.5,0,0),nrow=3)
M

#multiply matrix M(3x3) with a vector n0(1X3)
n<-M%*%n0
n

#multiply matrix M(3x3) with a vector n(3X1)
n1<-M%*%t(n)
n1

###############################################
# generate -
# a function that returns a matrix
# of iteratively added columns to a specified 
# vector n0. Each column is obtained as iterative 
# product of matrix M and vector n0. Returned  
# matrix has ngen+1 columns, where ngen is the  
# number of iterations that can be specified 
# by the user.
# input:      M -matrix
#             n0 - vector
#             ngen- number of iterations
# output:     a matrix of ngen+1 columns
###############################################
generate<-function(M,n0,ngen) {
#set up initial values for n and nrec
# to a vector n0
n<-nrec<-n0

#do in for loop 
#for each ngen
for (i in 1:ngen) {
#update vector n such that n is
#a product of the matrix  M(3x3) and 
#vector n(3X1) from previos iteration
n<-M%*%n
#update nrec such that current nrec is 
#a column binded from nrec from previos iteration
#and vector n from current iteration.
#(for the first iteration nrec is just a vector)
nrec<-cbind(nrec,n) 
}
#function outputs matrix nrec
return(nrec)
}
###############################################

#call the generate function 
#specify ngen to 10 iterations
ngen <- 10
nrec<-generate(M,n0,ngen=ngen)


# apply function
# runs the "sum" function
# to each of columns of the nrec matrix
# the second argument 2 refers to columns
#(if you want to sum up the row, the second 
#argument of the apply function should be 1)
(sums=apply(nrec,2,"sum"))

#plot columnwise sums of the nrec matrix
plot(0:ngen,sums,type="l")
#draw a line for the values from the first 
#row of the nrec matrix
lines(0:ngen,nrec[1,],lty=2)
#draw a line for the columnwise sum of the  
#first two rows from the nrec matrix
lines(0:ngen,apply(nrec[1:2,],2,"sum"),lty=3)


ev<-eigen(M)
nrecev<-solve(ev$vectors,nrec)
plot(0:ngen,nrecev[3,],type="l",lty=2)
plot(0:ngen,nrecev[1,],type="l",lty=2)
plot(0:ngen,nrecev[2,],type="l",lty=2)